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Abstract 



1 Introduction 



\^J Let us denote by D = K(x\, . . . , x n , d\, . . . , d n ) the Weyl algebra in n variables, 

J^ i where K is Q or C and di is the differential operator standing for Xi . We denote 

by D' = K(x m+ i, . . . , x n , d m+ i, . . . , d n ) the Weyl algebra in n — m variables, 
^^ . where m < n and D' is a subring of D. 

C*~) ' Let I be a holonomic left £>-ideal (|S]). The integration ideal of / with 

l/~j . respect to X\, . . . , x m is defined by the left £>'-ideal 

o, 

(I + diD + • • • + d m D) n D'. 

Oaku ([5]) gave an algorithm computing the integration ideal. This algorithm 
is called the integration algorithm for D-modules. The Grobner basis method 
p\ ' in D is used in this algorithm. 

i3 . We give a new algorithm computing not only generators of the integration 

ideal J but also Pq e / and Pi, ... , P m G D such as 

P = P a + d 1 P 1 + --- + d m P m 

for any generator P E J. Our algorithm is based on Oaku's one. We call these 
Pi, ... , P m inhomogeneous parts of P. As an important application of our algo- 
rithm, we can obtain inhomogeneous differential equations for definite integrals 
with parameters by using generators of the integration ideal and inhomogeneous 
parts. 
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For example, we compute an inhomogeneous differential equation for the 
integral A(x2) — J e~ Xl ~ XzX ^dxi. This is the case of m = l,n = 2. The 
annihilating ideal of the integrand f(xi, X2) = e~ Xl ~ X2Xl in D is / = (di + 1 + 
3x2x\, 82 + xf). The integration ideal of I with respect to x\ is J = (27:c2<9! + 
54x2^2 +6a; 2 + l) = (.P). The operator Pi = — (di+3<9i+3) is an inhomogeneous 
part of P. We apply the operator P to the integral A(x2) and obtain 
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3, 



P-A(x 2 )= d^Px ■ e- Xl - X2X i)d Xl 



P 1 -e 



-xi—x 2 x 1 



xi- b 



J ;ri=a 



(9x^xf - 3x 2 x{ - 6x 2 x 1 + l)e~ Xl - 

In this way, we get an inhomogeneous differential equation for the integral A{x2)- 
We will give an algorithm to compute inhomogeneous parts of the integra- 
tion ideal and give some examples. Other algorithms to compute differential 
equations for definite integrals are the Almkvist-Zeilberger algorithm (pQ, [TO] . 
[2]), the Chyzak algorithm ([4]) and the Oaku-Shiraki-Takayama algorithm ([7]). 
A comparison with these algorithms are also given. 

We implement our algorithms on the computer algebra system Risa/Asir 
I'll!]). They are in the program package nk_restriction.rr (|14). Packages 
Mgf un in Maple and HolonomicFunctions in Mathematica offers an analogous 
functionality, and are based on the Chyzak algorithm ([H], [13]). 

2 Review of the integration algorithm for D- 
modules 

We will review the integration algorithm for D-modules. We define the ring 
isomorphism J- : D — s> D satisfying 

J -a * (i<i<m) , m) = h ^ m \ . 

I Xi (to < 1 < n) \Oi (to < 1 < n) 

This map is called the Fourier transformation in D. 

The integration ideal of a left holonomic D-ideal / with respect to x\, . . . , x m 
is defined by the left D'-ideal J = (I + d x D + ■■■+ d m D) n D'. 

Algorithm 1 (Integration algorithm for D-modules, [6j, [§]) 

Input: Generators of a holonomic left D-ideal / and 

a weight vector w — (w\, . . . , w m , w m +i, . . . , w n ) such that wi, . . . , iu m > 
0,w m+ i = ••■ = iu„ = 0. 



Output: Generators of the integration ideal of / with respect to Xi 



■ • ) ^m- 



1. Compute the restriction module of the left Z?-ideal J- {I) with respect to 
the weight vector w. The details of the computation are as follows. 



(a) Compute the Grobner basis of the left D-ideal F{I) with respect 
to the monomial order <(- WjW ). Let the Grobner basis be G = 
{hi,..., hi}. 

(b) Compute the generic 6-function b(s) of J- {I) with respect to the 
weight vector (— w, w). 

(c) If b(s) has a non-negative integer root, then we set 
so = (the maximal non- negative integer roots). 
Otherwise, the integration ideal is and finish. 

(d) mi = ord(_ w , w )(hi), 

B d = {dl 1 — a*-* I hwi + ■■■ + i m w m <d} {d G N) , 

r = #{(ii,...,i m ) I hwi H h i m w m < s } = #S So . 

i 
(c) B = \J{hif, := dPfii | 0" G S S0 _ ro J, 
i=i 

S = {ft-i/3 := /lj/3 |a;i = ...=x m =0 I ^i/3 € B}. 

Here, h ll3 = J2d<*eB S0 9ad a {g a G £>')• 

2. Let (D') r be the left free L>'-modulc with the base J"" 1 (B ;So ), i.e. (D') r = 
^2d a eB D'x a . Regard elements in F~ l (B) as elements in the left D'- 
module (D') r . In other words, J r ~ 1 (hi^) = J2d a eB 9aX a (g a G £>') is 
regarded as an element in {D') r . Let M be the left D'-submodule in (D') r 
generated by J r ~ 1 {B). 

3. Compute the Grobner basis G of M with respect to a POT term order 
such that the position corresponds to a; = 1 is the minimum position. 
Output G' = GOD'. This set G' generates the integration ideal of /. 

We consider the following definite integral of a holonomic function f(xi, . . . , x n ). 

f(xi,...,x n )dxi ■•■dx m , R = JJ[ai,6»] 

1 i=i 

Let I — Anne/ := {P • / = | P G D} be the annihilating ideal of the 
integrand, and J be the integration ideal of /. For every p G J, there exist 
Pi , . . . , p m G D such that 



P-^ %>j G / 



r-l 



and we have 



p- A(x m+ i,...,x„) = / p-fdxi ■■■dx m = / y^(9iPi) ■ /rfxi ■■■rfx r , 



m „ 

V" / <%(Pi • /)dxi • • • dx m . (1) 

~1 Jfl 



Therefore, if we take an integration domain such that the right hand side of ([T]) 
equals to zero, we can regard the integration ideal as a system of homogeneous 
differential equations for the integral A(x m +ij ■ ■ ■ ; x n)- If the right hand side is 
not zero, the equation ([1} gives an inhomogeneous differential equations for the 
function A. 

3 Computing inhomogeneous parts of the inte- 
gration ideal 

In this section, we give a new algorithm of computing inhomogeneous differential 
equations for definite integrals. For the purpose, we must find an explicit form 
Pi (1 < i < m) in the equation ((T|) in the section [2j 

Theorem 1 Let J C D' be the integration ideal of a holonomic left D-ideal 
I . For any p G J, there exists an algorithm to compute differential operators 
Pi G D (1 < i < m) such that 

m 
p-J2 d iPi^ L ( 2 ) 

i=l 

Proof. We will present an algorithm of obtaining operators p t . By applying 
Algorithm [TJ we obtain a generating set {g\, . . . , gt} of the integration ideal of 
/. It is sufficient to compute inhomogeneous parts for each generator gj. From 
the step 3 of Algorithm [TJ gj can be expressed as gj — "^ J qjipJ-~ 1 {hip) where 
qji/3 G D. Then these qjip G D can be computed by referring the history of the 
Grobner basis computation in the step 3. Therefore, we have 

1 9 Yl Q^F' 1 ^) = 9j - \9j - Yl Qffl^&ip)) 

= 9j - Y q oifi ^ F ~ 1 (h i p) - F~ l (hip 

= 9i ~ Y^ 1 ! 3 \^~ 1 Q>'ip\x 1 =-=x m =o) - J" _1 (^/j) 

= 9j — 2^i 9j^^~ (hip\ Xl =—=x m =o — hip). 

Since each term of hip\ xl =...= Xm =o — hip can be divided from the left by either 
of xi, . . . , x m , each term of J- (hip\ xl =... =Xm= o — hip) can be divided from the 
left by either of d\ , . . . , d m . Thus we can rewrite 

m 

2 , qjipF^ (hip\ xl= ...= Xm =o — h.^) — 2_^ diPij. 

t=i 

Let us present our algorithm. 
Algorithm 2 



Input: Generators of a holonomic left ideal I C D and 

a weight vector w — (w± , . . . , w m , w m +i, . . . , w n ) such that W\, . . . , w m > 
0, w m +i = ■ ■ ■ = w n = 0. 

Output: Generators {gi, . . . , g t } of the integration ideal of / w.r.t. xi, . . . , x m and 
operators p^ € D satisfying gj — ^Z i=1 diPij € / for each generator gj 

a<j<t). 

1. Apply Algorithm [TJ 

2. Compute g^i^ satisfying gj = ^qjtp^ 1 ^^) by referring the history of 
the Grobner basis computation in the step 3 of Algorithm [TJ 

m 

3. Rewrite Rj := gj — ^2 qjipJ 7 ^ 1 ^^) to the form of Rj — \^ diPij- 



Output pij. 



i=l 



Example 1 [Incomplete Gauss's hypergeoemtric integral] 
We set 

F{x)= f t h - 1 (l-t) c - b - 1 (l-xt)- a dt. 
Jp 

We will compute a differential equation for the integral F(x). A holonomic ideal 
annihilating the integrand f(x,t) = t b_1 (l — t) c_6_1 (l — xt)~ a is 

I f =((-x 2 + x)d 2 x + ((-i + l)d t + (-a-b- l)x + c - l)d x - ab, 
(-t + l)xd x + {t 2 - t)d t + (-C+ 2)t + b - 1, (tx - l)d x + at) 

which is obtained by using Oaku's algorithm to compute the annihilating ideal 
of a power of polynomials. The generic b- function of 3^(1/) with respect to the 
weight vector w = (1, 0) (i.e. i's weight is 1 and a;'s weight is 0) is s(s — a + c— 1). 
We assume that a — c+ 1 is not a non-negative integer. Then the maximal non- 
negative integer root so of b(s) is 0. Therefore, the integration ideal of If with 
respect to t is 

((-x 2 + x)d 2 x + ((-a -b-l)x + c)d x - ab) = (P). 

The differential equation P • g = is Gauss's hypergeometric equation. The 
inhomogcneous part of P is dt(—t + l)d x . We apply P to the integral F(x) and 
obtain the inhomogeneous differential equation 



P- J f(x,t)dt=f {d t {t-l)d x )-f{x,t)dt 



df V 

(t-l)-£(x,t) 



We present the output for this problem by the program nk_restriction.rr 
on the computer algebra system Risa/Asir ([H]). We use the command 
nkjrestriction. integrationjLdeal to compute the integration ideal. The 
option inhomo=l make the system compute inhomogeneous parts and the op- 
tion param = [a,b,c] means that parameters are a,b,c. The sec shows the 
exhausting time of each steps. This example and next example are executed on 
a Linux machine with Intel Xeon X5570 (2.93GHz) and 48 GB memory. 



[1743] load("nk_restriction.rr") ; 

[1944] I_f=[-dx~2*x~2+(-dx*a-dx*b+dx"2-dx)*x-dx*dt*t-b*a+dx*c+dx*dt-dx, 

(-dx*t+dx)*x+dt*t~2+(-c-dt+2)*t+b-l,dx*t*x+a*t-dx] ; 

[(-x~2+x)*dx~2+((-t+l)*dt+(-a-b-l)*x+c-l)*dx-b*a, 

(-t+1) *x*dx+ (t"2-t ) *dt+(-c+2) *t+b-l , (t*x-l) *dx+a*t] 

[1945] nk_restriction. integration_ideal(I_i , [t,x] , [dt.dx] , [1,0] |param= 

[a,b, c] ,inhomo=l) ; 

— nd_weyl_gr :0.004sec(0.000623sec) 

— weyl_minipoly_by_elim :0sec(0.000947sec) 

— generic_bfct_and_gr : . 004sec(0.001922sec) 
generic bf ct : [ [1 , 1] , [s , 1] , [s-a+c-1 , 1] ] 

SO : 

B_{S0} length : 1 

— fctr(BF) + base :0sec(0.000277sec) 

— integration_ideal_internal :0sec(0. 000499sec) 
[[(-x~2+x)*dx~2+((-a-b-l)*x+c)*dx-b*a] , [ [[ [dt , (t-l)*dx]] ,1]]] 

Example 2 [F(x) = /"e"'- 1 ' 3 *] 

We consider the integral F(x) = L e~' _:E * dt. A holonomic ideal anni- 
hilating the integrand f(t,x) = e^*"** 3 is If = (d t + 1 + 3xt 2 ,d x + t 3 }. The 
integration ideal of If with respect to t is J = (27a; 3 9^ + 54x 2 d x + 6x + 1) = (P) . 
The inhomogeneous part of P is —dt(df + idt + 3). We apply P to the integral 
F{x) and obtain 

r-oo />oo 

P ■ / e-t-^dt = {d t {d 2 t + 3d t + 3)) • e- f - xt3 dt 

Jo Jo 

(d 2 + 3d t + 3) • e-*- x 

(-6xt + (1 + 3xt 2 ) 2 - 3 - 9xt 2 + 3)e" 

[1946] load("nk_restriction.rr") ; 

[2146] I_i=[dt+l+3*x*t~2, dx+t~3] ; 

[dt+3*t ~2*x+l , dx+t "3] 

[2147] nk_restriction. integration_ideal(I_f , [t,x] , [dt.dx] , [1,0] I 

inhomo=l) ; 

— nd_weyl_gr :0sec(0. 000526sec) 

— weyl_minipoly :0sec(0. 0002439sec) 

— generic_bf ct_and_gr : 0sec(0.001016sec) 
generic bfct : [ [1 , 1] , [s , 1] ] 

SO : 

B_{S0} length : 1 

— fctr(BF) + base :0sec + gc : 0.008sec(0.00691sec) 

— integration_ideal_internal :0sec(0 . 0003109sec) 
[[27*x~3*dx~2+54*x~2*dx+6*x+l] , [ [[ [dt ,-dt~2-3*dt-3] ] ,1]]] 

Theorem 2 We consider the following multiple integral, 

F(x m+1 ,...,x„) = / •••/ f(xi,...,x n )dxi---dx m (m<n). (3) 

J oi Jam 



Let I be a holonomic left D-ideal annihilating the integrand f(x\, . . . , x n ). There 
exists an algorithm to compute inhomogeneous differential equations for the mul- 
tiple integral F{x m -\-i, ■ ■ ■ ,% n ) from the holonomic ideal I. The algorithm is 
described below. 

For simplicity, we will explain the algorithm in the case of m — 2. We set 
F(x 3 ,...,x n ) = / f{x 1 ,...,x n )dx 1 dx 2 (2 < n), 

J a\ J a.2 

and will compute an inhomogeneous differential equation of F. 

Let J be a holonomic left D-ideal annihilating the integrand f(xi, ...,x n ). 
We compute the integration ideal J of J with respect to xi,x 2l i.e. 

J = (I + d x D + d 2 D) n D' (£>' = K{x 3 , ...,x n ,d 3 ,..., d n )). 

We take an element P E J. There exist P$ E I and Pi,P2 £ D such that 
P = P + d\P\ + d 2 P2 G D'. We apply the operator P to the integral F, and 
obtain 

PF= f \Prf\ Xl = bl -Prf\ Xl = ai )dx 2 + f ' (P 2 -f\ X2=b2 -P 2 -f\ X2=a2 )d Xl . (4) 

Let F\ , F 2 be the first term and the second term of the right hand side and let 
/ii fi be the integrand of F\,F%, 

To obtain a holonomic ideal annihilating the integral F\ , we must compute 
a holonomic ideal I\ annihilating the integrand f\. When the integrand f\ is 
the power of polynomial, we can use Oaku's algorithm to obtain the holonomic 
ideal I\ ([6]). In general case, we can compute the holonomic ideal I± from / by 
the following method. 

The ideal quotient I : Pi is holonomic and annihilates the function Pi- f. To 
obtain a holonomic ideal J\ annihilating P\ ■ f\ Xl= b 1 , we compute the restriction 
ideal oi I : P\ with respect to x\ —b\. Applying the same procedure for x\ = a% 
instead of x± = b\, we obtain a holonomic ideal J 2 annihilating P\ -f\ Xl=ai - Since 
J\ n J 2 is holonomic and annihilates /i(= P\ ■ /|xi=6i — P\ ■ /Ui=ai), we obtain 
J\ fl J 2 as I\ . 

We compute the integration ideal K\ of I\ with respect to x 2 , i.e. 

K x = (h + d^i) n D' (£>i = K(x 2 , x 3 , . . . , x n , d 2 , d 3 ,..., d n )). 

We take an element P (1) £ K x . There exist P (1) e h and P 2 (1) £ D x such that 
pW = P Q (1) + d 2 P 2 {l) . We apply F (1) to the integral F u and obtain 

P W . Fl = P 2 (1) ■ /iU 2=b2 - P 2 (1) • /iU 2=Q2 . (5) 

Applying the same procedure for I 2 instead of I\ , we can compute the annihilat- 
ing ideal I 2 of the integrand f 2 and the integration ideal K 2 of I 2 with respect 
to X\. 



By (JU and (JSJl, we obtain 

p(i) .p.p = p(i) . p 1 + p(i) . p 2; 

and can compute the first term of the right hand side. To compute the second 
term PW • F%, we compute K% : P^ 1 -* and take an element P^ in this ideal. 
Since pWpW g if 2 , we can compute p( 2 ')p( 1 ') • P 2 - Finally, we can obtain an 
inhomogcncous differential equation 

p(2)p(l)p . p _ p(2)p(l) . p i + p(2)p(l) . p 2 

Remark 1 Let 

£1 ■ F = gi, ■ ■ ■ , £ p ■ F — g p (1,6 P', g^ is a holonomic function ) 

be a system of inhomogeneous differential equations. When (£1, . . . ,£ p ) gener- 
ates the left holonomic ideal in D' , we call the system inhomogeneous holonomic. 
When 171 = 1, the output of the algorithm in Theorem [2] is inhomogeneous holo- 
nomic. Although the algorithm outputs a lot of inhomogeneous differential 
equations when P runs over the ideal J, it is an open question whether the 
output of the algorithm is inhomogeneous holonomic when m > 1. However, 
since the Oaku-Shiraki-Takayama algorithm gives holonomic output (see [7], 
section f4.3[) . we can obtain inhomogeneous holonomic differential equations by 
the following algorithm. 

Algorithm 3 

Input: Generators of a holonomic left ideal annihilating f{x\, . . . , x n ). 
Output: Generators of an inhomogeneous holonomic system for ([3]). 

1. Apply the algorithm in Theorem [5J 

2. Apply the Oaku-Shiraki-Takayama algorithm if the system obtained in 
the step 1 is not inhomogeneous holonomic. 

3. Merge the outputs of the step 1 and the step 2. 

4 Comparison of our algorithm with other algo- 
rithms 

4.1 The Almkvist-Zeilberger algorithm 

The Almkvist-Zeilberger algorithm (AZ algorithm, [TJ, [TU], [2]) is very fast, 
but works for hyperexponential functions. Our algorithm works for holonomic 
functions. The AZ algorithm is based on the method of undetermined coeffi- 
cients and Gosper's algorithm, and our algorithm is based on the Grobner basis 
method in D. 



4.2 The Chyzak algorithm 

The Chyzak algorithm ([3], @], [5]) is based on the method of undetermined 
coefficients and the Grobner basis method in the Ore algebra. By using the 
Ore algebra, the Chyzak algorithm can compute various summations and inte- 
grals like summations of holonomic sequences, integrals of holonomic functions 
and its ^-analogues. For the ring of differential operators with rational func- 
tion coefficients K(x)(d), the Chyzak algorithm is a generalization of the AZ 
algorithm and works for holonomic functions. The algorithm is often faster 
than our algorithm. But, when the algorithm returns higher order differen- 
tial equations or the number of variables are many, our algorithm is some- 
times faster. Here, we show only one example. We present these examples at 
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Example 3 [F(x, y) = J a xt+ y +tl0 dt] 

We set 

,6 1 

F(x,y)= / rrrdt. 

v ,y; J a xt + y + t w 

We will compute differential equations for the integral F(x,y). The following 
output is computed by our algorithm. It takes about 1.3 seconds. 

[2345] load("nk_restriction.rr") ; 

[2545] F=x*t+y+t~10$ 

[2546] Ann=ann(F)$ /* annihilating ideal of F~s */ 

0.052sec(0.0485sec) 

[2547] Id=map(subst, Ann, s, -1)$ /* substitute s=-l in Ann */ 

0sec(4.411e-05sec) 

[1569] nk_restriction. integration_ideal(Id, [t,x,y] , [dt,dx,dy] , [1,0,0] 

I inhomo=l) ; 

— nd_weyl_gr :0.012sec + gc : 0.008001sec(0.02009sec) 

— weyl_minipoly :0sec(0 . 001189sec) 

— generic_bfct_and_gr :0.016sec + gc : . 008001sec(0 . 02358sec) 
generic bf ct : [ [1 , 1] , [s , 1] , [s-9 , 1] ] 

SO : 9 

B_{S0} length : 10 

— fctr(BF) + base :0.044sec + gc : 0.024sec(0. 0674sec) 

— integration_ideal_internal :0.8321sec + gc : 0. 236sec(l .071sec) 
[ [9*x*dx+10*y*dy+9 , -10*dx~9-x*dy~9 , -9*dx~10+y*dy~10+9*dy~9] , 
[[[[dt,-t]] ,1] , [[[dt,-dy~8]] ,1] , [[[dt,-t*dy-9]] ,1]]] 

0.9081sec + gc : 0.28sec(l . 19sec) 

The following output is computed by the Chyzak algorithm (package Mgf un 
[12] ) on Maplel2. It takes about 50 seconds. 

with (Mgf un) : 

f :=l/(x*t+y+t~10) : 

ts :=time() : 

creative_telescoping(f , [x: :dif f , y : :dif f ] , t : :dif f ) : 

time()-ts ; 

49 . 583 



These computational experiments are executed on a Linux machine with 
Intel Xeon5450 (3.00GHz) and 32 GB memory. 

4.3 The Oaku-Shiraki-Takayama algorithm 

Although our algorithm gives inhomogeneous differential equations for defi- 
nite integrals, the Oaku-Shiraki-Takayama algorithm (OST algorithm, [7]) is 
for computing homogeneous differential equations annihilating definite integrals 
by using the Heaviside function and the integration algorithm. Since outputs 
are different, they are different methods. However, in most examples, outputs 
of our algorithm can be easily transformed to homogeneous systems. Thus, it 
will be worth making comparison between our method and the OST method. 

Let u(t, x) be a smooth function defined on an open neighborhood of [a, b] x U 
where U is an open set of R™ -1 . The Heaviside function Y(t) defined by Y(t) = 
(t < 0),Y(t) = 1 (t > 0). Then we can regard the integral of u(t, x) over [a, b] 
as that of Y(t — a)Y(b — t)u(t, x) over (— oo, oo), and the following holds. 

b />oo 

u(t,x)dt= Y(t-a)Y(b-t)u(t,x)dt 



Thus we can apply Algorithm [T] to obtain homogeneous differential equations. 
The paper [7] proposes the two methods 

(a) Method of using properties of the Heaviside function 

(b) Method of using tensor product in D-module 

to obtain differential equations annihilating the integrand of the right hand 
side. In the former case, the computation finishes without a heavy part because 
the procedure is only multiplication of polynomials. However, it is not known 
whether the output is holonomic. In the latter case, when an input is holonomic, 
an output is also holonomic. However, the computation is often heavy. We call 
the former OST algorithm (a) and the latter OST algorithm (b) in this paper. 
See Chap 5] for details. 

Let us show a relation of the outputs of OST algorithm and our algorithm. 
We consider v(x) = J °° e^ _t +t ^ x dt. OST algorithm (a) or (b) return the fol- 
lowing ideal 

( - 27x 3 <9 3 - 5Ax 2 d 2 + (Ax 3 + 3x)d x + Ax 2 - 3, 
27x 2 di + 135a;<9 3 + {-Ax 2 + 105)9^ _ \Q x d x - 8). 

On the other hand, Algorithm [2] returns the following ideal generated by P and 
its inhomogeneous part Q: 

(P) = {-27x 2 dl - Tlxd x + Ax 2 + 3), 
Q = d t (-9txd x + (-6t 2 + A)x + 3t). 
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This output yields 



p. v (x)= (-9txd x + {-6t 2 + 4)x + M) -e 



(-t 3 +t)x 



-Ax. 



t=o 



(6) 



Since the annihilating ideal of — 4x is (xd x — l,<9j), °P era tors (xd x — 1)P and 
d x P annihilate v(x). Although results of these algorithms are not coincide in 
general, these operators coincide outputs of OST algorithm (a) and (b) in this 
case. However, it seems that it is difficult to compute the right hand side of ^ 
from the output of OST algorithm. Moreover, in our algorithm we have only to 
do substitution process to compute for the integrals which has same integrand 
and another integration domain because our algorithm docs not depend on the 
integration domain. 

Tabic [1] shows the computing time of each part of Algorithm [5] and OST 
algorithm (a) , (b) . The entries with parentheses for inputs Vk mean that results 
for Vk were reused. For a comparison we show the computing time of Algorithm 
[TJ The experiments were done on a Linux machine with Intel Xeon X5570 
(2.93GHz) and 48 GB memory. 
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0.12 


0.0039 


V2 


0.15 


0.019 


0.17 


0.25 


5.10 


0.16 


5.26 


0.075 


V3 


19.91 


0.45 


20.36 


96.14 


24.54 


95.24 


119.8 


13.58 


V4 


26724 


28.33 


26752 


> 1 day 


1726 


> 1 day 




24003 


Vl 


(0.0042) 


0.0015 


0.0057 


0.0071 


0.47 


0.0050 


0.48 


n/a 


V2 


(0.15) 


0.027 


0.18 


1.56 


18230 


1.19 


18231 


n/a 


V3 


(19.91) 


1.62 


21.53 


3769 


848 


2802 


3650 


n/a 


Vi 


(26724) 


294 


27018 


> 1 day 


16231 


> 1 day 




n/a 



Vk(x) 



Uk(t,x)dt,Vk(x) 



,i / k 

/ Uk(t, x)dtwheieuk(t, x) = exp I — tx\ \(t —i ) 



Table 1: The comparison of the computing time (seconds) 

From the viewpoint of the computational efficiency, the computation time 
of Algorithm [5] increases more than that of Algorithm [1] for computing inhomo- 
geneous parts. That of OST algorithm increases because the input data of the 
integration algorithm becomes bigger differential operators by procedure (a) or 
(b). It seems that Algorithm [2] is faster than OST algorithm, since the compu- 
tation of inhomogeneous parts can be done by multiplication and summation of 
differential operators. However, to obtain homogeneous equation corresponding 
to OST algorithm output, we must compute annihilating ideals of inhomoge- 
neous parts. 
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